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Abstract 

The channeling of the ion recoiling after a collision with a WIMP in direct dark matter crystalline 
detectors produces a larger scintillation or ionization signal than otherwise expected. Channeling 
is a directional effect which depends on the velocity distribution of WIMPs in the dark halo of 
our Galaxy and could lead to a daily modulation of the signal. Here we compute upper bounds to 
the expected amplitude of daily modulation due to channeling using channeling fractions that we 
obtained with analytic models in prior work. After developing the general formalism, we examine 
the possibility of finding a daily modulation due to channeling in the data already collected by 
the DAMA/Nal and DAMA/LIBRA experiments. We find that even the largest daily modulation 
amplitudes (of the order of 10% in some instances) would not be observable for WIMPs in the 
standard halo in the 13 years of data taken by the DAMA collaboration. For these to be observable 
the DAMA total rate should be 1/40 of what it is or the total DAMA exposure should be 40 times 
larger. The daily modulation due to channeling will be difficult to measure in future experiments. 
We find it could be observed for light WIMPs in solid Ne, assuming no background. 
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I. INTRODUCTION 



The channeling effect in crystals refers to the orientation dependence of charged ion pen- 
etration in crystals. Channeling occurs when ions propagating in a crystal along symmetry 
axes and planes suffer a series of small-angle scatterings that maintain them in the open 
"channels" in between the rows or planes of lattice atoms and thus penetrate much further 
into the crystal than in other directions and loose all their energy into electrons. In dark 
matter crystalline detectors, a channeled ion recoiling after a collision with a WIMP (Weakly 
Interacting Massive Particle) would give all its energy to electrons, thus the quenching factor 
is Q ~ 1 instead of the usual Q < 1 for a non-channeled ion. Thus channeling increases 
the ionization or scintillation signal expected from a WIMP. The potential importance of 
the channeling effect for direct dark matter detection was first pointed out for Nal (Tl) 
by Drobyshevski [l{] and subsequently by the DAMA collaboration [2j in 2007. In 2008, 
Avignone, Creswick, and Nussinov [3] suggested that a daily modulation due to channeling 
could occur in Nal crystals, which would be a background free dark matter signature. Such 
a modulation of the rate due to channeling is expected to occur at some level because the 
'WIMP wind" arrives to Earth on average from a particular direction fixed to the Galaxy. 
Assuming that the dark matter halo is on average at rest with respect to the Galaxy, this 
is the direction towards which the Earth moves with respect to the Galaxy. Earth's daily 
rotation naturally changes the direction of the 'WIMP wind" with respect to the crys- 
tal axes, thus changing the amount of recoiling ions that are channeled vs non-channeled. 
This amounts to a daily modulation of the dark matter signal detectable via scintillation or 
ionization. 

Using analytic models of channeling which started to be developed in the 1960 's, shortly 
after the effect was discovered, mostly by Lindhard 4| and collaborators, we recently com- 
puted channeling probabilities as function of the recoil energy and initial direction q 



of a recoiling ion in different materials 
probability theory (see Eq. 5.13 in Ref. 



, |6] . We used a recursion of the addition rule in 
to find the probability x(Er, 4) that a recoiling 
ion enters into any channel in terms of the channeling fractions for single channels Xi{Er, q) 
that we computed (where the index i runs over all channels, both axial and planar). The 
channeling fractions for axial and planar channels are given in Eqs. 5.2 and 5.4 of Ref. jlj], 
respectively. 
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In our previous papers [5|, |6|, we also obtained the "geometric" channeling fraction 
-Pgeometric(-Ej?) in the crystals we studied, by averaging the channeling probability x(Er, q) 
over the initial recoil directions q (assuming an isotropic distribution in q) 



This integral was computed using the Hierarchical Equal Area iso-Latitude Pixelization 
(HEALPix) [7] of the recoil direction sphere (see Appendix B of Ref. jsj]). Here "geometric" 
refers to assuming that the distribution of recoil directions is isotropic. In reality, in a 
dark matter direct detection experiment, the distribution of recoil directions depends on the 
momentum distribution of the incoming WIMPs (see Section II). 

Fig. [TJa andCQb reproduced from Ref. show respectively upper bounds to some chan- 
neling fractions for single channels Xi{ERiQ) f° r Na recoils (with c = 1) and geometric 
channeling fraction of Na and I recoiling ions in a Nal crystal at room temperature for 
1 keV < Er < 20 keV. The parameter c mentioned in the figures is a number that we 
expect to be between 1 and 2, which regulates the importance of temperature corrections 
(for details see Re, fl). The channehng .actions are typically smaller for larger values of 
c thus setting c = 0, which is an unrealistic value, we get the largest upper bound to the 
channeling fractions that our calculations provide. In the figures we used c = and c = 1. 
Notice also that the results in the figures do not take into account dechanneling effects which 
should also decrease the channeling fractions (we do not know how to properly take into 
account these effects with our analytic methods). 

In this paper, we use the (upper bounds to the) channeling probability x(Er, q) and the 
actual differential recoil spectrum to compute the event rate, taking into account channeled 
and non-channeled recoils (see Section III, in particular Eqs. [T7] and [TS] and compare them 
with Eq. CD). We then use this rate to compute upper bounds to the amplitude of the daily 
modulation due to channeling expected in Nal crystals. In Section IV, we examine the 
possibility that such a daily modulation might be observable in the data accumulated by 
the DAMA collaboration. 




(1) 



3 



Na recoils, c=l 



0.05 




0.0015 



0.0000 L 



0.001 
5 X 10" 4 



1x10 



0.01 
0.005 



-4 



j i d. i i i i i i_ 

1. 1.5 2. 3. 5. 7. 10. 15.20. 
E (keV) 




5 10 15 20 



E (keV) 



FIG. 1: (Color online) Upper bounds to the (a) channeling fractions for single channels Xj(i?, q) 
of Na recoils for axial (black lines) and planar (green/gray lines) channels with c = 1, and (b) 
geometric channeling fraction -P ge ometric {E) of Na (solid lines) and I recoils (dashed lines) as a 
function of the recoil energy E for T = 293 K in with c = (green/gray) and c = 1 (black), always 
without including dechanneling. 

II. ANGULAR DISTRIBUTION OF RECOIL DIRECTIONS DUE TO WIMPS 

Consider the WIMP-nucleus elastic collision for a WIMP of mass m and a nucleus of 
mass M. The 3-dimensional "Radon transform" of the WIMP velocity distribution can be 
used to define the differential recoil spectrum as function of the recoil momentum q {^J 



where En is the recoil energy, dVL q = d<pd cos 6 denotes an infinitesimal solid angle around the 
recoil direction q = q/g, q = |q| is the magnitude of the recoil momentum, fi = mMj (m+M) 
is the reduced WIMP-nucleus mass, q/2fi = v q is the minimum velocity a WIMP must have 
to impart a recoil momentum q to the nucleus, or equivalently to deposit a recoil energy 
En = q 2 /2M, p is the dark matter density in the solar neighborhood, <To is the total scattering 
cross section of the WIMP with a (fictitious) point-like nucleus, and S(q) is the nuclear form 
factor normalized to 1. 

We concentrate here on WIMPs with spin-independent interactions, for which do is usu- 
ally written in terms of the WIMP-proton cross section a v [s| 



where /i p = mm p / {m + m p ) is the WIMP-proton reduced mass and A is the atomic number 




(2) 



cto = ^A 2 a, 



(3) 
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of the nucleus. We use the Helm form factor [1( 

,2 /'3ji(g#i, 



where 



S(q) = \F S i(q)\' 



qR 



smx 



cosx 



(4) 



(5) 



I' x 

is the first kind spherical Bessel function, R\ is an effective nuclear radius, and s is the 
nuclear skin thickness. Following Duda, Kemper, and Gondolo we set 



Ri 



c 2 + g 71-2(32 ~ 5s 2 , 



(6) 



and take s ~ 0.9 fm, a ~ 0.52 fm, and c ~ (1.23A 1 / 3 — 0.6) fm. These parameters have been 
chosen to match the numerical integration of the Two-Parameter Fermi model of nuclear 
density [11]. 



The Maxwellian WIMP velocity distribution with respect to the Galaxy, with dispersion 



a v and truncated at the escape speed t> esc is given by js] 



/wimp(v) 



for |v + Viab| < v esc , and zero otherwise, where 



No. 



erf 



cxp 



2v P . 



(v + V 



lab J 



2^ 2 



(7) 



exp 



V2aJ V 7T a v ^ [ 2(j2j- 
Here we are assuming the detector has a velocity Vi a b with respect to the Galaxy (thus 
— Vi a b is the average velocity of the WIMPs with respect to the detector). Vi a b is defined in 
terms of the galactic rotation velocity VcaiRot at the position of the Sun (or Local Standard 
of Rest (LSR) velocity), Sun's peculiar velocity V So i ar in the LSR, Earth's translational 
velocity VEarthRev with respect to the Sun, and the velocity of Earth's rotation around itself 
V Ea rthRot (see Appendix B), 



V 



lab 



V 



GalRot 



+ Vgolar + V 



EarthRev 



+ V 



EarthRot- 



(9) 



In this paper we take VcaiRot either 220 km/s or 280 km/s, as reasonable low and high values 
(as done in Ref 121]). which correspond to Vj a b either 228.4 km/s or 288.3 km/s, respectively 
(see Appendix B for details). Ref. 13] gives 100 km/s as the smallest estimate for the ID 
velocity dispersion, which corresponds to a 3D dispersion y/3 times larger, i.e. a v = 173 
km/s. Thus here we take a v either 173 km/s or 300 km/s j8j. 



In order to visualize the arrival directions of WIMPs, we will plot /wimp( v , v q ), the 
number of WIMPs per solid angle in the direction v in several figures. If we limit ourselves 
to the WIMPs with speed higher than v q , then 



/WIMP(V, V q ) 



/wimp(v)v dv. 



(10) 



The upper limit of the integral in Eq. [10] is such that |v + Vi a b| = f eS c and depends on the 
direction v, since (v + Vi ab ) 2 = v 2 + 2v v.V lab + VJ 



lab' 



-v.V lab + A /(v.V lab )2-V 1 2 ab + ^ s , 



and 



/wimp(v, v q 
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(11) 
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This integral can be solved analytically and the result is in terms of error functions, 



/wiMp(v, V qJ 



exp 



V 2 

Mab 

2a?, 



T 



^ [(v.V lab ) 2 + a 2 ] exp 



v.V 



lab J 



2a 2 



iV esc (27ra 2 ) 3 / 2 

erf ( V-Vlab + *Vax(v) \ _ erf / V-Vlab + ^ \ 

. V V2a v ) V V2a v J 

"max (v)(2v.Vi ab + Wmax(v)) 



+ (2(7, 



"max 

(v)) exp 



2^ 2 



+ (-V.Vi a b + V q ) exp 



t;g(2v.Vi ab + v q 
2al 



(13) 



The maximum of /wimp( v > v q) happens when v.Vi ab = — V^ ab , i.e. in the direction of the 
'WIMP wind" average velocity — Vi a b- Dividing /wimp(v, v q ) by this maximum we obtain 
a re-scaled distribution, a dimensionless number between and 1, which we plot in Fig. [2] 
(see the color scale/grayscale in the figure) on the sphere of velocity directions v using the 
HEALPix pixelization [7| (see also Appendix B of Ref. js|) for all WIMPs, which amounts 
to taking v q = 0. We took V[ ab = 288.3 km/s, and a v = 300 km/s or a v = 173 km/s for 
Fig. [2Ja or b respectively. 

For a truncated Maxwellian WIMP velocity distribution with respect to the Galaxy, 
truncated at the escape speed t> esc , the Radon-transform is jg] 



cxp 



[(q/2fi) + q.V 



lab 



iV esc (27ra 2 )V2 
if (q/2[i) + q.Vi ab < f eS c, and zero otherwise. 



2al 



— exp 



2^ 2 j 



(14) 




FIG. 2: (Color online) WIMPs number density per solid angle /wimp(v, v q ) (in Eq. [T3j) for all 
WIMPs (namely v q = 0) re-scaled to be a number between (black) and 1 (white) plotted on 
the sphere of velocity directions v using the HEALPix pixelization for Vj a b = 288.3 km/s and (a) 
a v = 300 km/s and (b) a v = 173 km/s. The arrow shows the direction of the average velocity 
of the WIMP wind, — Vi a b- The North and South celestial poles are also indicated. The color 
scale/grayscale shown in the horizontal bar between black and white corresponds to values between 
and 1 in increments of 0.05. 



The presence of q.Vi a b means that in order to compute the differential rate we need to 
orient the nuclear recoil direction q with respect to Vi a b. 



The maximum of /i a b(^"> q) m Eq. [H] happens when q.V 



lab 



-q/2fM, if v q = q/2fj,< V lah 



(or in the direction of — V lab otherwise). Thus, we can re-scale /i a b to obtain a dimensionless 
number between and 1, 



/) 



re— scaled 
lab 



exp 



[(g/2/x) + q.V lab ] 
2al 



exp 



2ol 




exp 



~^cs 

2^ 2 A 



(15) 



In Figs. [3] and H] we present side by side the WIMPs velocity distribution, for WIMPs 
which can generate a signal of a certain energy E, namely with speed above v q (left panels) 
and the Radon transform (right panels) of the recoils of energy E that WIMP collisions 
produce. 

In Fig. [3Ja and b we respectively plot /wimp(v, v q ) on the sphere of WIMP velocity 




FIG. 3: (Color online) (a) /wiMp(v,f g ) (in Eg. [T3j) re-scaled to be between and 1 plotted on 
the sphere of velocity directions v and (b) /i a b (re-scaled as in Eq. fl~5j) plotted on the sphere of 
recoil directions using the HEALPix pixelization for I recoils with Er = 10 keV, m = 30 GeV 
(thus v q = 304.6 km/s), V[ a b = 288.3 km/s and a v = 300 km/s. The arrow shows the direction 
of the average velocity of the WIMP wind, — Vi a b- The North and South celestial poles are also 
indicated. The color scale/grayscale shown in the horizontal bar corresponds to values between 
(black) and 1 (white) in intervals of 0.05. 

directions v and /i a b on the sphere of recoil directions (both re-scaled to be a number 
between and 1) using the HEALPix pixelization for I recoils assuming Vi a b = 288.3 
km/s, Er = 10 keV, a v = 300 km/s and m = 30 GeV. Fig. Hla and b show the same 
two distributions but for Na recoils and assuming a v = 173 km/s and m = 60 GeV (other 
parameters are the same). The color scale/grayscale plotted on the spheres indicate different 
values of the rescaled distributions: between (black) and 1 (white) in intervals of 0.05. 
In Fig. [3] the minimum WIMP speed required is v q = 304.6 km/s (I recoils), and since 
v q > Viab, the maximum value of /i^b~ scaLed , be. the maximum recoil rate, is in the direction 
of the "WIMP wind" average velocity, — V[ a b, which is shown with an arrow. In Fig. H] 
instead, v q = 196.7 km/s (Na recoils) and the maximum value of / la e b ~ scaled occurs when 
— q-Vjab = v q , i.e. when q is at an angle of 47° of — Vi a b- 
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FIG. 4: (Color online) Same as Fig. [3] but for Na recoils and assuming m = 60 GeV (so v q = 196.7 
km/s) and a v = 173 km/s (and all other parameters the same). 

III. DIFFERENTIAL ENERGY SPECTRUM 

Let p(E, Er, q)dE be the probability that an energy E is measured when a nucleus recoils 
in the direction q with initial energy Er, normalized so that 

J p(E,E R ,q)dE=l. (16) 

With our analytic approach we cannot estimate the importance of dechanneling mecha- 
nisms, such as the presence of lattice imperfections, impurities or dopants. Thus we disregard 
dechanneling, and assume that a recoiling nucleus can only either be channeled, in which 
case the measured energy is the whole initial recoil energy E = Er (first term in the fol- 
lowing equation) or not channeled, in which case the measured energy is E = QEr (second 
term) , 

p(E, E R , q) = X (E R , q)5(E - E R ) + [1 - X (E R , q)]6(E - QE R ). (17) 

The first term accounts for the channeled (unquenched) events and the second term for the 
unchanneled (quenched) events, and Q is the quenching factor. 
Using Eq. [T7] the differential energy spectrum, 

§=Ji£k mERA)dn < dE * (18) 
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can be written as 



dR 
dE 



dR 
dE 



+ 



dR 



dEftdVtq 



1 



dR 



E E ^ ^ //< ^' Q dE R dQ q 



e r =e/q 



dQ r , 



dR 



dE R dVt„ 



A . 1 

Eb=e ' *3 dE R dVt q 



e r =e/q 



dQ q , (19) 



where the differential recoil spectrum with subindex "U", which stands for "Usual" (i.e. 
when channeling is not taken into account) is 



dR 
dE 



1 dR 



Q dE R dVt q 



1 dR 



--E/Q 



Q dE 



R 



--E/Q 



(20) 



Defining q = y/2EM and using Eq. [2j the measured differential rate becomes, 



dR 
dE 



dR 
dE 



+ 



po"o 



S(d/VQ) 



S(q) I X (£,q)/iab(^,q)^ g 



q 



q 



q y — '- lab Vw^' 4 r"1' (21) 

Inserting Co from Eq. [3] in the above equation with the usual value for the mean local halo 
density p = 0.3 GeV/cm 3 , we can write the spin-independent detection rate of WIMPs in 
general for a crystal that may contain more than one element 



dR 
dE 



dR 
dE 



+ 1.306 x 10" 



events 



<7 44 



kg-day-keV Airm/ip 



S(q) / Xn{EA)h 



lab 



q 

2p r 



,q 



S(q/VQ~n) 



Xn(E/Q n ,q)f^b 



, q dVL q 



where 044 is the WIMP-proton cross section in units of 10~ 44 cm 2 , p, p and m are in GeV 
and J fiabdQq is in (km/s)" 1 . The sum is over the nuclear species n in a crystal, and 
Cn-, Xni Qn and \x n are the mass fraction, the channeling probability, the quenching factor 
and the reduced WIMP-nucleus mass for the element n, respectively. For example, for 
Nal crystals, as used in the DAMA experiment, we have C^a = ^Na/(^Na + Mi) and 
Ci = Mj/(Mj$ a + Mi), where M Na and Mi are the atomic masses of Sodium and Iodine 
respectively. 

The integrals in Eq. [25] cannot be computed analytically. We integrate numerically by 
performing a Riemann sum once the sphere of directions has been divided using HEALPix 7| 
(see also Appendix B of Ref. jf>]). HEALPix provides a convenient way of dividing the surface 
of a sphere into equal area sectors, and in our papers 0, 0] we use it for the first time to 
compute integrals over directions. 
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With the same notation, the usual rate is 



dR 
dE 
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events 



■x- 



cr 44 



kg-day-keV 47rm/i: 



/] Cn A\ 



Qn 
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lab 
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q dVt q 



(23) 



IV. DAILY MODULATION IN NAI CRYSTALS 



We present here the daily modulation amplitude due to channeling expected in Nal 
crystals for several WIMP masses and Na or I recoil energies. We assume that WIMPs have 
a truncated Maxwellian velocity distribution as in Eq. [7] with t> esc = 650 km/s. We use the 
upper bounds to channeling fractions for single channels Xi{Er, q) given in Ref. We take 
T = 293 K, the temperature of the DAMA experiment. 

The spin-independent detection rate of WIMPs given in Eq. [22] has a time dependence 
through the Radon transform f^. Notice that /i a b (see Eq. [14]) changes during a day through 
the (q.Vi a b) factor appearing in the exponent and the dependence of V lab on V EarthRot 
(see Eq. [9]). The expression showing the time dependence of q.Vi ab is given in Eq. IB13I 
(in Appendix B). During a day, VEarthRev which is responsible for the annual modulation 
changes too. Thus the rate does not return to exactly the same value after one day. For 
the cases we present in this paper, this difference is less than 10% of the total modulation 
amplitude in a day, and we did not correct for this effect. 



A. Relative Modulation Amplitudes 

Here we show the signal rate as function of time during a particular arbitrary Solar day 
(September 25, 2010). We define the relative signal modulation amplitude A s (taking into 
account the signal only) in terms of the maximum and minimum daily signal rate R s as 

3 = ~r 7p ' \ ' 

Ji'S-max "T -its— min 

The total relative modulation amplitude At is defined in terms of the maximum -Rr-max 
and minimum i?T-min total daily rates as 

A T = ^-max-i^-min (25) 
—max + R T —min 
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The total rate consists of signal plus background, Rt = R s + Rb- Assuming that there is no 
daily modulation in the background, -Rr-max — Rr-mm = R s -max — R s -mm, and At is related 
to A s as 



where the average total rate due to signal and background is Rt = (Rr-max + -Rr-min)/2 
and the average rate due to the signal alone is R s = (i? s _ max + R s -mm)/2- 

Exploring the parameter space of WIMP mass and WIMP-proton cross section for dif- 
ferent recoil energies we find that the relative modulation amplitudes A s can be large, even 
more than 10% for some combination of parameters. We explored the range of WIMP masses 
from a few GeV to hundreds of GeV for recoil energies between 2 keV and a few MeV. We 
show some examples in Fig. [51 where we plot the signal rate (in events/kg/day/keVee) as 
function of the Universal Time (UT) during 24 hours. We find that the largest A s happen 
when the signal is only due to channeling. This happens when there are no WIMPs in the 
galactic halo with large enough kinetic energy to provide the observed energy if the recoil 
is not channeled. The observed energies for which the rate is only due to channeling de- 
pend on the quenching factors Q, which are not well known. The smaller values of Q make 
channeling more important so we take = 0.2 [16] for Na and the usual Q\ = 0.09 for I. 




B. Statistical Significance 

The detectability of a particular amplitude of daily modulation depends on the exposure 
and background of a particular experiment. The former DAMA/Nal and the DAMA/LIBRA 
experiments (which we refer collectively as the DAMA experiment) have a very large cu- 
mulative exposure, 1.17 ton x year. However even with this large exposure, we find that 
the daily modulations we predict are not observable. To observe the daily modulation, the 
total number of events Nt {N s signal plus Nj, background events) over the duration of the 
experiment should be divided into two bins, the "high-rate" bin with A^_ max events and the 
"low-rate" bin with iV T _ min events, so that N T = A^_ max +iVV_ mm . For the daily modulation 
to be observable at, say, the 3a level one should have 



A T = A S (R S /R T ), 



(26) 




N T . 



—max 



— Nt 



—mm 



A T N T > 3a ~ 3y/N T /2, 



(27) 
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FIG. 5: Signal rate (in events/kg-day-keVee) as function of the Universal Time (UT) during 24 
hours for m = 10 GeV, 12 GeV and 15 GeV for different energies. The parameters used are 
a v = 300 km/s, (^Na = 0.2, Qi = 0.09, o p = 2 x 10 _40 cm 2 , c = 1 for temperature effects, a crystal 
temperature of T = 293 K and Vi a b = 228.4 km/s (top row) or 288.3 km/s (bottom row). 



where a 2 ~ N T /2 because, with a small modulation, on average iV" T „ max ~ N T _ min ~ N T /2. 
In principle there are other errors associated with identifying the "high-rate" and "low-rate" 
bins which we do not include here. Thus we are underestimating the errors. 

If the detector exposure is MT in kg-day and we take bins of width AE in keVee, 
then iV T _ max = R T _ max MTAE/2, N T _ min = R T _ min MTAE/2, N T = R T MT AE and 
N s = R S MTAE, where the rates are in events/kg-day-keVee. Thus (N s /Nt) = (R s /Rt) 



and using Eq.EH A T = A S (N S /N T ). Thus the condition in Eq. M becomes A S N S > 3^/N T /2 
which implies 

N^/N T > 9/(2A 2 s ), (28) 



or 



Rj/R T > 9/(2A 2 s MT AE). (29) 
The total rate of the DAM A experiment at low energies 4 keVee < E < 10 keVee is 



Rt — 1 events/kg/day/keVee 



3- 



This rate is much larger than the signal rates we predict 



and is, therefore, dominated by background. With this value of Rt, Eq. |2"91 becomes 

Q 

s s 2MT AE kg day keVee' 1 ' 
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We choose here a bin AE ~ 1 keVee, narrow enough to assume the signal rate to be 
constant in it and compatible with the energy resolution of DAMA. The energy resolution 
of DAMA is c E (E) = (0.448 keVee)^£/keVee + (0.0091)E ~ 1 keVee at low energies (l5|. 
We consider the significance of the highest signal-to-noise energy bin that we found through 
inspection. With the cumulative exposure of DAMA, the condition in Eq. [30] for relative 
daily modulation amplitudes A s observable at 3cr is 

R s A s > 3.2 x 10~ 3 events/kg/day/keVee, (31) 

or 

-Rs-max — -Rs-min > 6.4 x 10~ 3 events/kg/day/keVee . (32) 

For observability at the na level we should multiply the right-hand side of Eq. [32] by (n/3). 
Even the largest relative daily modulations we find, shown in Fig. [5] are not observable in 
the DAMA data according to Eq. I3"2l 

The examples which we show here are for small WIMP masses and recoil energies. For 
large masses the value of a p must be chosen in the region of the cross section and mass plane 
where XENON10/100 and CDMS impose a p to be smaller by four orders of magnitude than 
for light WIMPs. This amounts to corresponding smaller signal rates and (-R s - ma x — R s -mm) 
differences. For small WIMP masses and large energies, v q is large and there are no WIMPs 
with the speed required for Na or I recoils. Thus, only small WIMP masses and recoil 
energies result in high modulation amplitudes. 

Fig.[5]shows the signal rate during 24 hours for three different WIMP masses m = 10 GeV, 
12 GeV and 15 GeV and different energies E between 2 and 15 keVee. The other relevant 
parameters are a v = 300 km/s, a p = 2 x 10 _40 cm 2 (close to the DAMA and CoGeNT 



regions |16l4l8l|). c = 1, T = 293 K and two values of V[ a b, 228.4 km/s (top row) and 288.3 
km/s (bottom row). Recent bounds, e.g. those from XENON100 jl9j ]. impose smaller values 
of a p . In any event, changes in a p are easy to take into account because A s is independent 
of a p and the rate is just proportional to it, R s ~ a p . 

We found the relative amplitude A s to be as large as 12% in the examples shown in Fig. El 
but even those large values are not observable according to Eq. [32] (even at the la level). 
With the choice of V[ a b = 228.4 km/s (top row of Fig. [5]) we get a signal rate difference 
-Rs-max — -Rs-min of 0.56 x 10~ 3 events/kg/day/keVee for m = 10 GeV and E = 10 keVee 
(in this case v q = 454.8 km/s and 790.5 km/s for channeled Na and I recoils, respectively), 
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3.17 x 10~ 4 events/kg/day/keVee for m = 12 GeV and E = 12 keVee (which corresponds to 
v q = 441.6 km/s and 732.9 km/s for Na and I channeled recoils, respectively), and 4.25 x 10 -4 
events/kg/day/keVee for m = 15 GeV and E = 15 keVee (for which u g = 430.6 km/s and 
670.6 km/s for Na and I channeled recoils, respectively). With the choice of = 288.3 
km/s (bottom row of Fig.E]), -R s - ma x~ -Rs-min is 0.77x 10 -3 events/kg/day/keVee for m = 10 
GeV and E = 10 keVee (one of the energies shown), 2.95 x 10~ 4 events/kg/day/keVee for 
m = 12 GeV and E — 12 keVee, and 0.58 x 10~ 5 events/kg/day/keVee for m — 15 GeV and 
-E 1 = 15 keVee. Because the minimum WIMP speeds v q are large in these examples, a smaller 
velocity dispersion of the WIMP distribution leads to smaller rates (since a smaller amount 
of WIMPs have velocities larger than v q ). So the signal rate difference -R s _ max — R s -min is 
even smaller for smaller values of a v . 

The left-bottom panel of Fig. [5] shows the signal rate as function of UT for m — 10 GeV 
and V| a b = 288.3 km/s for several energies between 2 keVee and 12 keVee. The rate decreases 
but A s increases with increasing energy and the best conditions for observability happen at 
some energy where neither the rate nor A s are very small. The rates for low energies between 
2 keVee and 6 keVee are dominated by the usual (i.e. non-channeled) rate and the daily 
modulation is due purely to the change in WIMP kinetic energy in the lab frame as the 
Earth rotates around itself. The rates for energies above 8 keVee (green/gray lines) are 
purely due to channeling, i.e. the usual rate is zero. For intermediate energies, 6 keVee to 8 
keVee, the usual and channeled rates both contribute and thus the daily modulation is due 
to both the channeling effect and the daily change in the usual rate. For E = 2, 4, 6, 8, 10 
and 12 keVee, the values of R s - max — R s -mm given in events/kg/day/keVee are respectively 
4.3 x 10~ 4 , 0.5 x 10" 3 , 0.92 x 10" 3 , 2.8 x 10~ 4 , 0.77 x 10~ 3 and 0.52 x 10" 3 . Notice that 
for all the energies shown the difference in rate is similar, but the largest A s values happen 
at energies above 8 keVee, for which the rate is only due to channeling. The channeling 
daily modulation amplitude increases as the ratio of the velocity dispersion to the average 
speed of the WIMPs that contribute to the signal (i.e. with v > v q ) decreases. This ratio is 
small and thus A s large for large values of v q . Notice that the phase of the modulation due 
to channeling depends on the orientation of the crystal with respect to the Galaxy and the 
phase of the modulation in the usual rate does not, which would allow to distinguish both 
effects, if they were observable. The case of m — 10 GeV and E — 6 keVee has the largest 
rate difference, but is not observable at 3<r according to Eq. |32] (not even at the ler level). 
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c=l, E=3.8 keVee, m=5 GeV, V GalRot =220 km/s 



c=0, E=3.8 keVee, m=5 GeV, V GalRot =220 km/s 
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FIG. 6: Signal rate as function of UT during 24 hours for E = 3.8 keVee and m = 5 GeV, with 
VLab = 228.4 km/s, <r„ = 300 km/s, <r p = 2 x 10- 40 cm 2 , and Q Na = 0.2, Q x = 0.09 for (a) c = 1 
and (b) c = 0. The daily modulation is not observable in both cases. 



Choosing a p = 4 x 10~ 40 cm 2 (still within the DAMA allowed region but not compatible with 
the recent XENON100 result) results in a rate difference of 1.84 x 10 -3 events/kg/day/keVee 
for this case which would not be observable even at the la level. 

Finally, we would like to compare our results with those obtained in Ref. Q by Creswick 
et al. They found a relative daily modulation amplitude A s =0.85% (their definition of 
amplitude differs by a factor of 2 from ours, so they quote 1.7%) for 5 GeV WIMP mass and 
3.8 keVee measured energy (in which case v q = 471.2 km/s and 936.6 km/s for channeled Na 
and I recoils, respectively. There are no WIMPs with the speed required for I recoils, thus 
only Na recoils are possible). In order to compare our calculation with theirs, we compute 
the signal event rate as function of time for c = 1, T = 293 K (temperature corrections are 
not included in the calculation of Creswick et al.) and choosing all the other parameters very 
close to those used in Ref. [3], i.e. V[ a b = 228.4 km/s and a v = 300 km/s. A WIMP mass 
of 5 GeV is outside the region of parameter space compatible with the annual modulation 
reported by DAMA [l^]. Since A s does not depend on a p , we choose an arbitrary value of 
a p = 2 x 10~ 40 cm 2 to plot the signal rate as a function of UT (the upper bound given by 
TEXONO and CoGeNT [20] is five times larger, ff p <lx 10~ 39 cm ). Our result is shown 
in Fig. El a. We find A s =0.16% {R s - m * x - R s -mm = 4.4 x 10" 6 events/kg/day /keVee). 
Even when we consider the extreme choice of c = to compute temperature effects (an 
unrealistic value for which the channeling fractions are larger) with the same parameters, 
we get A s = 0.14%. This case is shown in Fig. [6jb. 



16 



C. Future Prospects for DAM A and other Experiments 



The daily modulation might be detectable in other experiments with smaller background 
or WIMP halo components with a smaller dispersion such as streams or a thick disk. The 
amplitude of the daily modulation increases as the WIMP velocity distribution is narrower 
i.e. for larger values of the average velocity and smaller values of the velocity dispersion of 
the detectable WIMPs (which is not a v ), i.e. those with velocity larger than v q . This is easy 
to understand since as the dispersion increases more channels are available for channeling 
of the recoiling ions. In the limit in which the velocity distribution would be isotropic with 
respect to the detector, the daily rotation would not introduce any difference in the rate 
due to channeling. Having a large relative signal modulation amplitude A s is not sufficient 
for observability. In Eq. [32] what is important is (A s R s ) = (i? s _ max — i? s _ mm )/2. However, 
the condition in Eq. [32] was derived considering the total rate in the DAMA experiment, 
which is dominated by background. For an experiment where the background is negligible, 
i.e. Rt = R s + Rb — R s , we can derive a different observability condition (at the 3a level) 
from Eq. [231 

R 8 A 2 S = A s (R s „ max - R s - min )/2 > 9/(2MT AE). (33) 

This condition might be easier to satisfy in future experiments. 

One could ask which is the maximum level of total rate with the current DAMA exposure 
that would be needed to make the signal daily modulation observable. Inserting the current 
exposure of DAMA (1.17 ton year) in Eq. [29] we have 

(A s Rsf/Rr > 1-05 x 1(T 5 events/kg/day/keVee, (34) 

which using A S R S = (R s -m^ ~ R s -mm)/2, becomes 

R < (Rs—m&x -Rs— min) (35) 

T 4.2 x 10- 5 events/kg/day /keVee' 1 ' 

Even in the case with the highest rate difference we found, i.e. -R s - m a X — -Rs-min = 0.98 x 1CT 3 
events/kg/day/keVee (the m = 10 GeV, E = 6 keVee, V| a b = 288.3 km/s example shown in 
the bottom-left panel of Fig. [5]) observability would require 

R T < 0.023 events/kg/day /keVee, (36) 

roughly 1/40 of what is now. 
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We could ask instead what exposure would be needed with the current total rate 
in the DAMA experiment to make the daily modulation observable. Setting — 1 
events/kg/day/keVee in Eq. [2H1 we obtain 

MTAE > 9 = 18 (37) 

(events/kg/day/keVee) 2 (A a R s ) 2 (i? s _ max - i? s _ min ) 2 ' 

Again, for the case with the highest rate difference we found (m = 10 GeV, E = 6 keVee 
and Vi a b = 288.3 km/s) and with AE ~ 1 keVee we would require an exposure 40 times 
larger, 

MT > 51.3 ton year. (38) 

We have computed the daily modulation due to channeling in other material such as 
Ge, solid Xe and solid Ne, and we find that it will be very difficult to observe. For light 
WIMPs the cross section can be larger than for heavier ones without violating experimental 



bounds, Or, = 10 39 cm 2 



'v 



201 ] and this favors the detection of the daily modulation. We 



find that for a WIMP mass m = 5 GeV the daily modulation due to channeling may be 
observable in solid Ne if the signal would be above threshold and assuming no background. 
The geometric channeling fraction reaches a maximum at around 10 keV for solid Ne Q], 
thus the largest modulation amplitude happens at that energy. For example for a solid Ne 
detector operating at 23 K at Gran Sasso, for E — 10 keV, assuming Q^ c = 0.25 [2 if . c = 1 
and with velocity distribution parameters o v = 300 km/s and V[ a b = 228.4 km/s we find 
R S A 2 S = 3.68 x 10~ 5 events/kg/day/keVee. Using Eq. [33] we find that the exposure needed 
to observe this modulation at 3<r is MT = 0.33 ton year. For the same parameters but for 
m = 7 GeV and o v = 2 x 10~ 40 cm 2 (parameters compatible with the possible dark matter 
signal found by CoGeNT and with DAMA according to Ref. [22J|), we find R S A 2 S = 7.2 x 10~ 7 
events/kg/day/keVee, and the exposure needed is MT = 17 A ton year. The usual rate is 
zero in both cases, and the modulation is just due to channeling. The signal rate during 
24 hours and the required exposures for the two cases are shown in Fig. [7] and Table [H 
respectively. 

We intend to further explore the observability of a daily modulation in future experiments 
for different halo models in future work. 
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c=l, E=10 keVee, m=5 GeV, o- p =10" 39 cm 2 c=l, E=10 keVee, m=7 GeV, o- p =2xl(T 40 cm 2 

5 10 15 20 04 5 10 15 20 
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FIG. 7: Signal rate as function of UT during 24 hours for a solid Ne detector operating at T = 23 
K at Gran Sasso for E = 10 keVee, Q = 0.25, c = 1, a v = 300 km/s, Vj a b = 228.4 km/s and for (a) 
m = 5 GeV and cr p = 10~ 39 cm 2 , and (b) m = 7 GeV and a v = 2 x 10~ 40 cm 2 . 
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0.0098 
0.0096 
0.0094 
0.0092 
0.0090 
0.0088 



TABLE I: Observability in solid Ne detector 



Case 


Op (cm 2 ) 


MT (ton year) 


m = 5 GeV 


10 -39 


0.33 


m = 7 GeV 


2 x 10- 40 
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Appendix A: Crystal Orientation 

We need to orient the crystal with respect to the laboratory. We define a reference frame 
fixed with the laboratory and orient its axes so that the xy plane is horizontal, the x-axis 
points North, the |/-axis points West, and the z-axis points to the zenith. We denote its 
unit coordinate vectors as Af, VV and Z, respectively. We also define the crystal frame with 
X, Y, Z cartesian axes fixed with the crystal. The unit coordinate vectors of the crystal 
frame are X, Y and Z. 

We now want to connect the laboratory frame to the crystal frame. Let the standard 
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orientation correspond to the configuration in which X = Af , Y = VV, and Z = Z. We start 
with the crystal in the standard orientation, and we turn it into any other orientation X, 
Y, Z. In this new orientation, each of the unit coordinate vectors of the crystal frame can 
be written in terms of unit coordinate vectors of the lab frame, 

X = a x Af + Px VV + 7x Z, 
Y = a Y Af + Py VV + 1y Z, 

Z = a z Af + p z VV + j z Z, (Al) 

where Pi and 7« are the "direction cosines" between the two sets of cartesian coordinates 
of the lab and crystal frames, for i = X,Y, Z. For example, the coordinate vector X of 
the crystal has a particular angle with each of the lab frame coordinate vectors Af, VV, Z. 
Let ax be the angle between X and Af, bx the angle between X and VV, and cx the angle 
between X and Z. The direction cosines of the unit vector X are given by, 

ax = cos ax = X • Af, 
Px = cos bx =X - VV, 

7x = cosc x = X-Z. (A2) 

We can find the direction cosines for Y and Z unit vectors in a similar way. From these 
definitions it follows that a.i aj + Pi Pj + ji jj = 5ij where i,j = X,Y, Z. We prefer using 
direction cosines over Euler angles because the direction cosines can easily be measured for 
any known orientation of a crystal in a laboratory, whereas it may be difficult to specify the 
Euler angles. 

Eq. IA1I gives the transformation from the lab frame to the crystal frame. We can also 
find the lab coordinate vectors in terms of the crystal coordinate vectors, 

Af = a x X + a Y Y + a z Z, 
VV = p x X + Py Y + p z Z, 

Z = lx X + 7 y Y + lz Z. (A3) 

In the results we show in this paper, we took ax = Py = lz = 1 and all the other oti, Pi 
and 7j equal to zero. Choosing a different orientation for the crystal does not change the 
average rate, but A s may change by a factor of 2 for Nal depending on the orientation of 
the crystal. The observability condition is still not satisfied. 
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1. Lab to equatorial transformation 



To connect the laboratory frame to the equatorial coordinate frame, we recall the defini- 
tion of the geocentric equatorial inertial (GEI) frame: its origin is at the center of the Earth, 
its x e -axis points in the direction of the vernal equinox, its y e -axis points to the point on 
the celestial equator with right ascension 90° (so that the cartesian frame is right-handed), 
and its z e -axis points to the north celestial pole. We denote its unit coordinate vectors as 
x e , y e , and z e . We want to find the transformation formulas from the laboratory frame to 
the GEI frame. 

This transformation can be achieved by two successive rotations. The first rotation is by 
an angle of (90° — Ai a b) counterclockwise about the laboratory y-axis to align the new x'y' 
plane with the plane of the celestial equator. Here Aj a b is the latitude of the laboratory in 
degrees, with northern latitudes taken as positive and southern latitudes taken as negative. 
With this rotation, the new z'-axis points to the north celestial pole. The second rotation is 
by an angle (15ti a b + 180) degrees clockwise about the new z'-axis to bring the x'-axis in the 
direction of the vernal equinox. Here ti a b is the laboratory Local Apparent Sidereal Time 
(LAST) in hours (the LAST is the hour angle of the vernal equinox at the location of the 
laboratory). One has 

tlab = ^GAST + WIS, (A4) 

where £qast is the Greenwich Apparent Sidereal Time (GAST) in hours and /i a b is the 
longitude in degrees measured positive in the eastward direction (e.g. Zi a b = +110° for 110° 
E and h ab = -110° for 110° W). 

The current local apparent sidereal time for any specified longitude Zi a b can be 
computed online, for example on the website of the US Naval Observatory at 



http: / /tycho. usno.navy.mil/] sidereal.html (accessed Sept 19, 2010). As an alternative, one 



can use the following formula 23|, |24j for the Greenwich mean sidereal time (which differs 



from the Greenwich apparent sidereal time by less than 1.2 seconds, completely negligible 
for our purposes), 

t GAST = (101.0308 + 36000.770 T + 15.04107 UT)/15, (A5) 

where 

IMJDJ - 55197.5 , A , 

T = ± J - . A6 

36525.0 v ; 
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Here UT is the Universal Time in hours, |_MJDJ is the integer part of the modified Julian 
date (MJD), which is the time measured in days from 00:00 UT on 17 November 1858 (Julian 
date 2400000.5). Note that T is the time in Julian centuries (36525 days) from 12:00 UT on 
1 January 2010 to the previous midnight. At 12:00 UT on 1 January 2010, the Julian date 
is 2455198, and the MJD is 55197.5. Also the the 15.04107/15 in Eq. \M\ corrects from solar 
time (UT) to sidereal time. Sidereal day is shorter than Solar day by 3.9 minutes. In this 
paper, all our results are computed for the particular arbitrary day of 25 September 2010, 
for which T = 0.00729637. 

Note also that UT is different from coordinated Universal Time (UTC) which is the time 
scale usually used for data recording. UTC is atomic time adjusted by an integral number 
of seconds to keep it within 0.6 s of UT. For our purposes the difference between UT and 
UTC is negligible. 

Taking into account the two rotations explained above, one can find the transformation 
equations of the unit vectors, 

x e = - cos(t lab ) sin(Ai ab )A/ r - cos(A lab )2 + sin(t° ab )VV, 
y e = - sin(t° ab ) sin(Ai ab )A r - cos(Ai a b)i" - cos(t lab )W, 
z e = cos(Ai ab )A/' + sin(Ai a b)2, (A7) 

where t° ab = 15ti a b is the laboratory LAST converted to degrees. 

As a check, for a laboratory on the equator at local sidereal time 0, i.e. Ai a b = 0° and 
£° ab = 0°, one has x e = Z, y e = — VV, and z e = Af; six sidereal hours later at the same 
laboratory, i.e. Ai a b = 0° and i° ab = 90°, one has x e = VV, y e = Z, and z e = Af; for a 
laboratory at the South Pole (Ai a b = —90°), using the direction of the Greenwich meridian 
in place of the "North" axis Af so that the local sidereal time at the South Pole by convention 
coincides with the Greenwich sidereal time, one has x e = Af, y e = — VV, and z e = — Z at 
£° ab = 0° and x e = VV, y e = Af, and z e = — Z at i° ab = 90° . All of these are correctly given 
by Eq. E3 

The formulas in Eq. IA7I can be inverted, and the transformation from the equatorial 
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FIG. 8: (Color online) Earth's sphere in the equatorial frame (x e ,y e ,z e ) specified with black arrows. 
The laboratory frame (N,W,Z) specified with blue/dark gray arrows is also shown. 

frame to the lab frame is achieved: 

J\f = - sin(Ai ab ) [cos(t° ab )x e + sin(t° ab )y e ] + cos(Ai ab )z e , 
W = sin (*°ab)x e - cos(*i ab ) y e , 

Z = cos(Ai a b) [cos(t° ab )x e + sin(t lab )y e ] + sin(Ai ab )z e . (A8) 

The latitude and longitude of Gran Sasso are Aj a b = 42.45° and /i a b = 13.7°, respectively. 

Fig. [8] shows the laboratory frame (Af, VV, Z) and the equatorial coordinate frame 
(x e ,y e ,z e ) plotted on the Earth's sphere at UT = using Eq. IA8I 

2. Equatorial to galactic transformation 

To connect the equatorial frame to the galactic coordinate frame, we recall the definition 
of the galactic coordinate system: its origin is at the position of the Sun, its x 9 -axis points 
towards the galactic center, its y 9 -axis points in the direction of the galactic rotation, and 
its points to the north galactic pole. 

For the epoch of January 1950.0 the transformation from the equatorial frame (x e , y e , z e ) 
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to the galactic frame (x 9 ,y 9 ,z 9 ) is given by [251 ]: 



x 9 = x e (-0.06699) + y e (-0.8728) + z e (-0.4835), 
y 9 = x e (0.4927) + y e (-0.4503) + z e (0.7446), 

z 9 = x e (-0.8676) + y e (-0.1883) + z e (0.4602). (A9) 

The transformation from the galactic frame to the equatorial frame is given by 

x e = ± g (-0.06699) + y g (0.4927) + z g (-0.8676), 

y e = x 9 (-0.8728) + y 9 (-0.4503) + z g (-0.1884), 

z e = x g (-0.4835) + y g (0.7446) + z g (0.4602). (A10) 

The change of Eqs. IA9I and IA10I from the epoch of January 1950.0 to 25 September 2010 is 
small and would not affect the final results in this paper. 

Appendix B: Laboratory motion 

The velocity of the lab with respect to the center of the Galaxy can be divided into four 
components (as in Eq. EJ): V Ga iRot, V So i ar , V Ea rthR ev and V Ea rthRot- 

We take V GalRot = 220 km/s or 280 km/s ^ Solar = 18 km/s V EaithRev = 29.8 
km/s and VEarthRot = (0.465102 km/s) cos Ai a b, where Ai a b is the latitude of the lab. Values of 
KjaiRot — 220 km/s or 280 km/s results in V[ a b = 228.4 km/s or 288.3 km/s, respectively (see 
Appendix B.5 for the equation of Vi a b). Thus, Vi a b is dominated by the galactic rotation 
velocity. 

We need to compute q • Viab, where q is given in the crystal reference frame (q = 
qx X + qy Y + qz Z). Therefore, we need to also write Vi ab in the crystal frame. We have, 

q ■ Vi a b = q ' VcalRot + q • Vsolar + q - VsarthRev + q ' VEarthRot- (Bl) 

We will compute each term on the right-hand side of Eq. IB1I individually. 
1. Galactic rotation 



The velocity of the galactic rotation VcalRot is defined in the galactic reference frame, 

VcalRot = VcalRotyg, (B2) 
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where VcaiRot is the galactic rotation speed (i.e. the local circular speed), and y g is in the 
direction of the galactic rotation. Following Ref. jjjj , we take VcaiRot = 220 km/s or 280 
km/s. Using the conversions in Eq. IA9[ we can write y g in the equatorial reference frame in 
terms of (x e ,y e ,z e ). Then, we use Eq. IA7l to transform from the equatorial frame to the lab 
frame (Jv, VV, Z), and finally we use Eq. IA3l to transform from the lab frame to the crystal 
frame (X,Y,Z). 

Thus, we can use Eq. IA3I to write VcaiRot in terms of the crystal frame coordinates, and 
compute q • V Ga iRot, 

q • VcaiRot = <?xVGalRot,X + ?yVGalRot,Y + 9Z^GalRot,Z- (B3) 

We have 

q ■ V GalRot = y Ga iRot| ( [-0.4927 cos^J + 0.4503 sin(t° ab )] sin(A lab ) + 0.7446 cos(A lab ) 

(a x q x + a Y q Y + a z q z ) + ( 0.4927 sin(^ ab ) + 0.4503 cos^) j (f3 x q x + (3 Y q Y + fa 
+ 

ilxQx + IyQy + Izqz) }■ (B4) 



[0.4927co S ( Cb ) -0,503 S i„(C b )lco S (A„ b ) +0.744.^) 



Eq. IB4l has a time dependence through t° ab and would be responsible for any daily modulation 
in the rate. 

2. Solar motion 

The velocity of the Sun's motion in the galactic rest frame is, 

Vgolar = Ux g + Vy g + Wz g , (B5) 



where (U, V,W)q = (11.1,12.2,7.3) km/s [26j. Using Eq. \A9\ we can transform from the 
galactic frame to the equatorial frame, and using Eq. IA7I we can transform from the equato- 
rial frame to the lab frame. Then we can use Eq. IA3I to write Vg i ar in terms of the crystal 
frame coordinates. 
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Thus, we can compute q • Vsoiar as 
q ■ Vsoiar = M (1.066 km/s) cos(£° ab ) + (16.56 km/s) sin(t° ab )] sin(Ai ab ) + (7.077 km/s) cos(Ai ab 
(a x qx + azyCtY + Q-zqz) + ( - (1.066 km/s) sin(t^ b ) + (16.56 km/s) cos(t° ab 
(/W + PyQy + Pzqz) + (- [(1-066 km/s) cos(t° ab ) + (16.56 km/s) sin(t° ab )] cos(A 
+ (7.077 km/s) sin(Ai ab ) J {^ x qx + Ivqv + Izqz)- (B6) 



lab 



Clearly, Eq. IB6I has a time dependence through i° ab and would be responsible of any daily 
modulation in the rate. 

3. Earth's revolution 

The velocity of the Earth's revolution around the sun is given in terms of the Sun ecliptic 



longitude A(t) as [28 1 



V E arthRev = Vs(X(t)) [COS /3(x) SHl(A(t) - X x )* g 

+ cos f3(y) sin(A(t) - A^)^ + cos 0(z) sin(A(t) - \ z )z g ], (B7) 

where = 29.8 km/s is the orbital speed of the Earth, V^(A(t)) = V e [l — esin(A(t) — A )], 
e = 0.016722, and Ao = 13° + 1° are the ellipticity of the Earth's orbit and the eclip- 
tic longitude of the orbit's minor axis, respectively, and = (—5°. 5303, 59°. 575, 29°. 812) 
and A; = (266°. 141, -13°.3485, 179°.3212) are the ecliptic latitudes and longitudes of the 
(x.g,y g ,Zg) axes, respectively. 



The Sun's ecliptic longitude \(t) can be expressed as (p. 77 of Ref. [27J] and Ref. [281 ]) 



A(t) = L + (1°.915 - 0°.0048T ) sin g + 0°.020 sin 2g, (B8) 

where L = 281°. 0298 + 36000°. 77T + 0°.04107£/T is the mean longitude of the Sun corrected 
for aberration, g = 357°. 9258 + 35999°. 05To + 0°.04107£/T is the mean anomaly (polar angle 
of orbit). 

Using Eq. \A9\ we can transform from the galactic frame to the equatorial frame, and 
using Eq. I A 71 we can transform from the equatorial frame to the lab frame (J\f, \V,Z). Then 
we can use Eq. IA3I to write Vsoiar in terms of the crystal frame coordinates. 
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Thus, we can compute q • VEarthRev as 
q • V EarthRev = V e (A(t)) j [ - cos(t° ab ) sin(Ai a b)^4(0 - sm(t° ab ) sin(Ai ab )B(i) + cos(Ai ab )C(t)] 
(a x qx + ayqv + a z qz) + [sin(£ lab )„4(t) - cos(t° ab )£(£)] ((5 x qx + fiyqv + PzQz) 
+ [cos(t° ab ) cos(Ai ab )^l(t) + sin(t° ab ) cos(Ai ab )#(£) + sin(Ai ab )C(£)] (j x qx + Ivqv + Izqz) 

(B9) 

where 

A{t) = (-0.06699) cos/3(x) sin(A(t) - \ x ) + (0.4927) cos 0(y) sin(A(t) - X y ) 

+ (-0.8676) cos p(z) sin(A(t) - A z ), 
B(t) = (-0.8728) cos/3(x) sin(A(t) - \ x ) + (-0.4503) cos/%) sin(A(t) - X y ) 

+ (-0.1883) cos/3(z) sin(A(t) - A 2 ), 
C(t) = (-0.4835) cos/3(x) sin(A(t) - A*) + (0.7446) cos 0(y) sin(A(t) - X y ) 

+ (0.4602) cos p(z) sin(A(t) - \ z ). (B10) 

Eq. IB9I has a time dependence through t° ab and A(t) and would be responsible for any 
daily modulation in the rate. 



4. Earth's rotation 

Finally, we want to compute VEarthRot, the velocity of Earth's rotation around itself. We 
have 

V E arthRot = -^RotEqCOsAiabVV, (Bll) 

where VptotEq is the Earth's rotation speed at the equator, and is defined as ViiotEq = 
2ttRq/(1 sidereal day). The Earth's equatorial radius is R 9 = 6378.137 km, and one sidereal 
day is 23.9344696 hr= 86164 s. therefore V RotEq = 0.465102 km/s. 

Using Eq. IA3I to write VV in terms of the crystal frame coordinates, we can easily find 
q ■ V EarthRot as 

q • V E a rt hRot = - VRotEq cos Ai ab (Pxqx + fivqy + fizqz) • (B12) 

There is no time dependence in Eq. IB12| because it is written in the crystal frame, and both 
the lab and the crystal are rotating with the Earth. 
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5. Total Velocity 



Now we can insert Eqs. IB41 IB61 IB9l and IB12l into Eq. IBll to compute q • Vi a b- Inserting 
the values of V® = 29.8 km/s, e = 23.439° and pRotEq = 0.465 km/s, we have (in km/s): 



q V 



lab 



cos(t° ab ) A(t)+sm(tl h ) B(t) 



sin Ai ab + C(t) cos A lab \ {a x q x + a Y q Y + ot z qz) 



+ <^ sin(^ ab ) A(t) + cos(t° ab ) B{t) - 0.465 cos A lab \ ((3 x q x + fi Y q Y + fi z q z 



+ 



cos(t° ab ) A{t) - sin(t° ab ) B{t) 



cos Ai ab + C{t) sin A lab \ {ixQx + IyQy + IzQz) , 

(B13) 



where 



A(t) = 0.4927 y GalRot - 1.066 km/s + {V @ (X(t))A(t), 
B{t) = 0.4503 V GalR ot + 16.56 km/s - (V e (X(t))B(t), 
C{t) = 0.7445 V GalRot + 7.077 km/s + (V e (A(*))C(t). 



(B14) 
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